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A generalized Kuramoto model of coupled phase oscillators with slowly varying coupling matrix 
is studied. The dynamics of the coupling coefficients is driven by the phase difference of pairs of 
oscillators in such a way that the coupling strengthens for synchronized oscillators and weakens 
for non-synchronized pairs. The system possesses a family of stable solutions corresponding to 
synchronized clusters of different sizes. A particular cluster can be formed by applying external 
driving at a given frequency to a group of oscillators. Once established, the synchronized state 
is robust against noise and small variations in natural frequencies. The phase differences between 
oscillators within the synchronized cluster can be used for information storage and retrieval. 
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I. INTRODUCTION 



The mechanisms of adaptation and learning in natural and artificial systems are the subject of paramount interest 
in neuroscience. Most of the experimental evidence point to the synaptic modification as the physiological basis of 
the long-term memory [0|2) . Hebb ||] was the first who suggested that the synaptic coupling between two neurons is 
enhanced if both neurons are simultaneously active. Recent experimental data suggest that the synaptic connections 
are strengthened (long-term potentiation) or weakened (long-term depression) depending on the precise relative timing 
between pre-synaptic and post-synaptic spikes [^,^||^]. 

Many current theoretical studies of learning in neural systems are based on the quasi-static Hopfield model |^ 
utilizing McCuUoch-Pitts neuronal units or simple "integrate-and-fire" models of neurons coupled via synapses which 
accumulate the pre-synaptic activity (spikes from other neurons) into a growing membrane potential of a neuron [Q. 
The latter generates its own spike ("fires") once the membrane potential exceeds a certain threshold. The synaptic 
strength is in fact not fixed, but changes in response to external stimuli and/or inter-neural dynamics. Various models 
based on this general picture have been shown to yield differential Hebbian learning rules. 

A complimentary approach to modeling the neuronal activity involves replacing neurons by periodic oscillators. 
Indeed, rhythmic activity plays important role in many neuronal systems and functions, including central pattern 
generators, visual and olfactory systems, etc. Neural networks can often be described as networks of coupled 

oscillators, so the role of relative spike timing is played by the phases of individual oscillators (see, for example, 
pd]-p^). Close to the Hopf bifurcation, the dynamics of the array of coupled oscillators can be described by the 
Landau-Stewart equations for complex amplitudes of the oscillators. Furthermore, one can reduce the dynamics of 
the oscillators to the pure phase dynamics by assuming that the coupling is weak, and the amplitude of oscillations 
is slaved to the phase variations. In fact, the phase dynamics description may still be valid even away from the 
Hopf bifurcation point, where the Landau-Stewart description is not applicable. This idea was first put forth by 
Winfree and later developed by Kuramoto ||l6| and others The original Kuramoto model, in which the 

coupling coefficients were chosen fixed and equal, allowed for the elegant analytical treatment within the mean field 
approximation. It was shown that the system of globally coupled oscillators with non-identical frequencies exhibits a 
second-order phase transition as the coupling strength is increased above some critical value. 

In subsequent work, more complex coupling functions were considered ||l^,|l^,|l9| . It was shown that depending on 
the choice of the connectivity matrix, synchronized states with non-trivial phase relationship among the oscillators 
can emerge. Networks of coupled phase oscillators with appropriately tuned coupling matrix were shown to have neu- 
rocomputing properties similar to those of classical Hopfield networks, with an important distinction that memorized 
patterns are stored in the form of relative phases of synchronized oscillators rather than in static equilibria. 

In those studies, the connectivity matrix, however complicated, was imposed externally to achieve the desired 
network dynamics. In fact, in order to achieve certain behavior (for example, to memorize an image), in a system 
of N neurons, N'^ coupling coefficients have to be specified. This would present significant difficulties in operating 
such a system at large N . In biological neural networks this task is not assigned to any central control unit, but is 
performed in a distributed manner via the mechanism of synaptic plasticity, i.e. long term potentiation or depression 
of synapses in response to the dynamics of pre-synaptic and post-synaptic neurons External stimuli directly 

affect only the dynamics of individual neurons and not their connections, while the synaptic plasticity is a result 
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of the inter-neuronal interactions. This evolution of connectivity matrix have been introduced in Ref. jT^ for the 
Landau-Stewart model of coupled oscillators, however no analysis of the joint dynamics of the network of oscillators 
and the connectivity matrix has been presented. In this paper, we study the simplest possible model of this process, 
based on the generalization of the Kuramoto model. In particular, we assume that the coupling coefficient for a link 
between two oscillators is a slow function of the phase difference between the two oscillators. 

The paper is organized as follows. In Section || we introduce the ge nera lized Kuramoto model with additional 



equations describing the slow evolution of the coupling matrix. In Section III, a simple case of two coupled oscillators 
is investigated. We show that for large enough separation of time scales of fast and slow motions, the system is 
bistable: depending on initial conditions, oscillators can be either synchronized or non-synchronized. In Sections 
^,[v|, we study the existence and stability of the family of synchronized cluster solutions in the continuum limit 



N ^ oo. In Section VI, we show that the generalized Kuramoto system can serve as a very simple model of learning 



and memory. Information is stored within the synchronized cluster in the form of phase differences between oscillators. 



Section VII summarizes our conclusions. 



II. MODEL 

Let us consider the dynamics of an ensemble of N coupled phase oscillators 

1 ^ 

•^i = - — ^ K,jF{(f)i - (f>j). (1) 

Here uji are natural frequencies and are phases of individual oscillators, F{(j)) is a 27r-periodic coupling function, 
and Kij is the N x N matrix of coupling coefficients. 

In order to include the mechanism of adaptation into model (0) , we assume that an element of the coupling matrix 
describing interaction between two oscillators, i and j, is controlled by the following equation, 

X,, =e(G(0, (2) 

where G{4>) is a 27r-periodic function of its argument. The particular form of the coupling function F{(j)) and adaptation 
function G{(p) can be derived from the underlying equations describing the dynamics of oscillators through the 
reduction to the phase description. Following Kuramoto [|l6|, we choose the simplest possible periodic function 
F(4)) — sin (/). Furthermore, we choose the adaptation function G{(j)) — a cos (/). This choice implies that the coupling 
coefficient grows fastest for two oscillators which are in-phase and decays fastest for out-of-phase oscillators. There 
are certain indication in the biological literature (see Ref. |^) that the maximal LTP and LTD correspond to small 
positive and negative phase shifts between pre-synaptic and post-synaptic neuronal oscillations, but we are going to 
ignore this complication in this model for the sake of simplicity. With the chosen functions, our model becomes 



0, = - ^ sin((/)j - (3) 
i=i 

kij = e(Q!Cos((/)j - (f>j) - Kij). (4) 



For small e, the dynamics of the coupling coefficients is slow, and one can expect that for two oscillators which are 
non-synchronized, the driving term acos{(j)i — 4>j) is oscillating around zero, and the resultant coupling coefficient 
is also oscillating around zero and is small 0(e). On the other hand, a pair of synchronized oscillators produces a 
constant non-zero driving for the corresponding coupling coefficient, and therefore for large enough a, the connection 
between the two will be strengthened. Thus, we anticipate a multi-stable behavior, when a group of oscillators can 
either be stably synchronized or non-synchronized depending on the initial conditions for their coupling coefficient. 
To illustrate this point, we first consider a simple case of two oscillators symmetrically coupled through a single link. 

III. TWO COUPLED OSCILLATORS 

In this case model reduces to a set of two coupled equations, 

(j) = Acj — K sin 

k = €{acos(j)- K). (5) 
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Here (f> = — (f>2, K = K12 = K21, and Aw = lji — 102- By introducing new variables t = l^ujt^K = K/lS.uj^a = 
a/Aw, e = e/Aoj, system (H) is simplified to 

(j) — \ — K svtv(f), 

K = e{a cos (f) — K) . (6) 

The nuUclines of this system K = acos(/) and K = l/sin(/> intersect at large enough a > etc — 2. In this case, 
within each period of the relative phase 0, there are four intersections corresponding to two stable and two unstable 
fixed points (see Figures H,b,c. For small e, the dynamics of the system can be separated into fast and slow motions. 
If initial value Kq > 1, the phase variable rapidly approaches the quasi-static value arcsini^'o"^ corresponding to 
the nuUcline, without any significant change of K. Then, depending on whether Kq is greater or smaller than the 
value Ku corresponding to the unstable fixed point, the solution cither slowly approaches the stable fixed point Kg 
(synchronized regime), or it reaches the minimum of the nuUcline and then branches off to infinity along the fast 
trajectory. The same occurs if Kq < 1. This corresponds to the regime of non-synchronized oscillations. For larger 
values of e, the dynamics can no longer be reduced to fast and slow motions. Moreover, for any fixed a > dc, 
there exists some critical value of ic such that at e > Ec there are no unbounded solutions corresponding to the 
non-synchronized motion. On the other hand, if a < cic, there are no fixed points corresponding to the synchronized 
solutions, and the oscillators are desynchronized for any values of e. Figure ||b,c shows the structure of the phase 
plane for system (||) for e > Ic and e < Cc, respectively. Figure |l| shows the dependence of ic vs. d. As expected, 
ic diverges as a ^ Q:c+. This line together with vertical line d = etc divides the parameter plane into three regions. 
For d < dc^ there can only be non-synchronized solutions (Fig.|^,a), for d > dc,e > Ic there are only synchronized 
solutions (Fig.^b), and for d > dc,i < f-c, there is a bistable state when both desynchronized and synchronized 
regimes are stable (FigJ|,c). 

The results of this section can be applied to a more generic case of many coupled oscillators. Indeed, as e and/or 
a is increased, the spontaneous clustering usually begins from synchronization of pairs of neighboring oscillators. For 
this process, the coupling between these oscillators and other oscillators can be neglected. Symbols in Fig. |l| show 
the results of numerical experiments with a population of 50 oscillators uniformly distributed within a range Afi = 1 
starting from small randomized initial coupling Kij . The bifurcation to random parings of oscillators which precedes 
the cluster formation, occurs close to the theoretical line predicted for a simple case of two coupled oscillators. 
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FIG. 1. Phase diagram ec{& for the model of two coupled oscillators (g). Symbols show the results of numerical simulations 
with 50 uniformly distributed oscillators. Circles correspond to the non-synchronized regime, and stars indicate an onset of 
synchronization of pairs of oscillators, which usually precedes spontaneous cluster formation. 
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IV. MANY COUPLED OSCILLATORS - STATIONARY STATES 



In this section we return to the dynamics of the orig inal model (|),(|)- As was pointed out in Section ||, for small 
e, in the asymptotic regime the coupling coefficients are either non-stationary but small (for two oscillators operating 
at different frequencies) or fixed for two oscillators which are synchronized. So, in the limit e — > 0, we can neglect all 
Kij for non-synchronized oscillators, and define A'y = Q;cos(0i — for synchronized oscillators. Substituting this 
into Eq.(^, we obtain for oscillators within a cluster, 

= ^» - ^ X! sin[2(0i - (f>j)]. (7) 

where unlike Ref. |]l6| summation only applies to oscillators within the same cluster, and 4'i — for oscillators outside 
synchronized cluster. It is easy to see that Eq. (|^) is formally equivalent to the Kuramoto model for the double phases 
2(j)i. 

One can investigate this model in the mean-field limit iV — > oo using the order parameter method [ p^ . In this case, 
the complex order parameter reads 

2iip 



re^''= hm N'^Y e''^^ = / 



n{iP)e''^d^. (8) 



where ip = cj) — Qt, 17 is the frequency of the synchronized cluster, n^ip) is the phase distribution of the oscillators in 
the population which belong to the cluster. The phases of the oscillators inside the cluster are given by 

^.^g+^arcsinf ^(^'~") V (9) 

and the phase distribution can be easily expressed via the frequency distribution of the synchronized oscillators 

nU\ = [ '^^^('^ -Q)^^) cos(2(^ - Q)) for B<2(^i,-B)<A 

^ ' I 0, outside ^ ' 



where A^B determine size and position of the synchronized cluster, (— 7r/2 < B < A < tt/2). Substituting (|I0|) into 
Eq.(||), we obtain the equation for the magnitude of the order parameter r 



/A 
3 sinx -f f7j cosx e"da; = 1. (11) 

The frequency spread inside the synchronized cluster is determined by formula 



A17c = — (sinA-sinS), (12) 

which follows directly from Eqs.(j|). Thus, this system exhibit degeneracy so that a variety of synchronized states 
can be formed depending on initial conditions. It is particularly evident for the uniform distribution of oscillator 
frequencies within frequency interval Afi, g{uj) = Af2~^. In this case, the cluster size drops out, the limits of 
integration are symmetric, B = —A, and we arrive at 



a 



2MI 



A 

cosx e"da; = 1, (13) 

-A 



which reduces to an algebraic equation for phase width of the cluster A^ 

A+isin(2A) = ^. (14) 
2 a 

The phase width of the cluster A depends only on a/Afi and does not depend on the cluster size r. This dependence 
is shown in Fig.|^. For large a, A « Ari/a. At a ^ ao = 47r~'^Ar2, A 7r/2. The frequency spread inside the 
synchronized cluster is equal to Afic = arsiwA (cf. (|l^). At a < ag, there can be no stationary synchronized 
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clusters. These findings are confirmed by direct numerical simulations of Eqs.(|^),(Q). In Fig.|^ the phase width of the 
cluster is shown for several values of a/AJl, and in Fig.^ the evolution of an unstable cluster at a = 1.25 < ao is 
shown. The parameters of simulations were N — 50, Ail — 1, e — 0.002, and the initial cluster width was chosen 
A^c 0.5. 
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FIG. 3. Phase spread of the cluster ^4 as a function of a/ An from Eq.Q. Dots indicate the results of numerical simulations 
of Eqs.(|i),(|) with 50 uniformly distributed oscillators with An = l,e = 0.002. 




5 10 15 20 25 30 35 40 45 50 
Oscillator index 

FIG. 4. Gray-scale plot of the evolution of instantaneous oscillator frequencies (pi for an unstable cluster at a = 1.25 
Parameters of the simulation are A'^ = 50, e = 0.002, Ail = 1. 



V. STABILITY OF THE CLUSTERED STATE 



As we have seen in the previous Section, the order parameter r which characterizes the cluster size, can take any 
value from to Tmax (for the uniform frequency distribution the latter corresponds to the whole population being 
synchronized, Tmax = Ail/ a sin A). However, this is true only for small enough e. Just as in the case of two coupled 
oscillators (Sec. HI), at sufficiently large e > Ec, the multistability disappears, and the only stable state at a > etc is 
the one involving the whole population of oscillators, AQc = Ail. In order to determine ec, we analyze the following 
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problem. Suppose we have a synchronized cluster oscillating at frequency and characterized by the order parameter 
T. Let us consider a single oscillator with frequency loq adjacent to the cluster (wq = + Ar2c/2), and study its 
dynamics neglecting its interaction with all other oscillators and the the feedback influence of the oscillator on the 
cluster. Clearly, since the frequency difference between this oscillator and the cluster is smaller than that for any 
other non-synchronized oscillator, the adjacent oscillator will be the first to get entrained at Cg. The equation for the 
relative phase ipo = (po — fit of this oscillator reads 



1 ^= 

-^X,sin(V'i - V'o), (15) 



i=l 

and the equation for the coupling coefficient between the oscillator and the i-th member of the cluster, 

ki ^ e{acoa{il;i - ^o) ~ Ki). (16) 

The summation in Eq.ffq) is carried only over synchronized oscillators. The phases ipi of the synchronized oscillators 
are determined by Eq.(P)^ 

In the continuum limit for the uniform frequency distribution, these equations become 

^^ = ^ + —j^^^ K{i^) cos(2V') sin(V' - Vo)rf^, (17) 
k{iP) = e{a cos{4> - V'o) - K{iP)). (18) 
Eq.(pT[) can be rewritten in the form 

A57 ar 

V-o = + ^(Pcos^o - QsinVo), (19) 



where 



.A/2 

P= K{ijj) sin ip cos 2'ipd'ip, (20) 

J-A/2 
.A/2 

Q= K{iJj)cos->pcos2ipdil). (21) 



Equations for P, Q can be obtained using Eq.dTl), 



P = e (aP_ sin ~ P) , (22) 
Q = e(aP+cosi/'o-Q), (23) 



whereP± = ^sin A±{jA+ ^sm2 A). The set of equations (19),(p2|),(p3[) exhibits the dynamics similar to Eqs.(|^) for 
two coupled oscillators. At a < ao — An^^Ail, the only attractor of the system is the solution with periodic P and 
Q and unbounded phase 4> which corresponds to "non-entrainment" of the oscillator by the cluster. For the phase 
variable taken modulo 2tt, this solution corresponds to a limit cycle encircling the phase cylinder. At larger a > ao, 
two fixed points appear, one of which (stable) corresponds to the entrainment of the oscillator by the cluster. At 
e > eo(a), the limit cycle disappears, and the entrained state is the only stable state of the system. Thus, the critical 
line 60(0;) separates the region of cluster stability with respect to entraining additional oscillators. This line is shown 
in Figure ^. This line should be compared with the critical line for pairing of two neighboring oscillators (Figj^). 
Note that we have to re-scale e and a back into original e and a using e = eN, a = 2a/ Ail. For large N, this line 
lies below eo{a), which indicates that the random pairing outside the cluster occurs before the main cluster loses its 
stability. This conclusion agrees with our numerical simulations. 
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FIG. 5. Parameter plane (a,e). A cluster is stable at a > qo and e < eo. For a < qo, clusters do not exist, for a < ao 
and e > eg, cluster size increases spontaneously. Dashed line shows the stability limit for synchronization of two neighboring 
oscillators, for a > ac and e > ec, random pairing of oscillators outside the main cluster occurs. 



VI. LEARNING AND MEMORY 



Information can be stored within the cluster in the form of relative phases of the oscillators in the synchronized state. 
Initial learning can be accomplished via driving a group of oscillators by external signals with identical frequency loq 
and magnitude Kq, but with different phases 



1 ^ 

X! ^"-1 ~ ^j) + ^0 sin(0i - ujot - cj)°), 



e{a cos( 



(24) 
(25) 



If Kq ^ 1, every oscillator in the group will be entrained to the external frequency uiq. So, the external driving signals 
force the oscillators to oscillate in synchrony and, via slow synaptic dynamics Eq. (|2^) , form a tightly coupled cluster. 
If the phases 0° are taken to be or tt, the relative phases — (f>j within the cluster will be also close to or tt, and 
so the coupling coefficients Kij will approach zLa. After the learning is completed after t ^ 0(e^^) time, the external 
signal is disconnected {Kq — > 0). If the external frequency luq = fl, the mean frequency of oscillators within the 
cluster, the cluster remains synchronized at the same frequency, moreover, the phase relations among the oscillators 
are effectively preserved. Even significant random fluctuations affecting the dynamics of individual oscillators, do not 
change the robust structure of the synchronized cluster, see Fig.^. 

For small e, after a cluster is formed, it becomes very stable with respect to external noise. Furthermore, it remains 
stable with respect to random variations of their natural frequencies. We let the cluster form similarly to described 
above, and then at some time ti, we perturbed the natural frequencies of all oscillators, uji = tUi + Auj^i, where S,i 
are i.i.d. Gaussian random variable with variance 1. We can quantify the robustness of the cluster by measuring the 
correlation between the synchronized cluster structure before and after the random fluctuations of natural frequencies 
have been applied. The correlation was calculated using the formula 



C 



(26) 



where Pi{qi) is the binary variable describing the state of the oscillator before(after) perturbing the natural frequencies: 
Piili) = 1 if the oscillator is entrained in the main cluster and Pi{qi) = otherwise. Symbols V and A denote Boolean 
addition and multiplication, respectively. It is easy to see that C — I when cluster remains unchanged, and C = if 
all oscillators leave the cluster. 
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FIG. 6. Relative phase distribution of oscillators within the cluster with respect to the "middle" oscillator i = 25. The 
cluster is initialized by external forcing of oscillators 15-37 at frequency uio — = 10.0 for < t < to — 20. At f > to, the 
external forcing was turned off, but the phase distribution within the cluster survived. At t > ti = 75, oscillators were driven 
by random external fluctuation of magnitude 0.01. Despite of this forcing, the phase pattern "memorized" by the cluster, was 
preserved. 

Figure ^ shows correlation C as a function of the perturbation magnitude Auj for several values of a. As can be 
seen in the Figure, the cluster becomes more robust with increase of a. 
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FIG. 7. Stability of cluster measured with the correlation of the cluster before and after perturbation of natural frequencies. 
The different curves are for several values of a, e increases the stability of the cluster only slightly and falls of parallel to the 
curves for similar a. 



VII. CONCLUSION 



In summary, we have demonstrated that a generalized Kuramoto model with additional equations describing slow 
dynamics of the coupling matrix can be used to elucidate the underlying mechanism of synaptic plasticity. The slow 
dynamics leads to the multi-stability: synchronized clusters of different sizes and with different phase relationships 
among oscillators can be stabilized. The phase differences among the oscillators can be used as a way of storing and 



retrieving the information in this system. One natural hmitation of the reduced phase description of the oscillators is 
that all oscillators are assumed to be in the excited ("firing") state. However, in natural systems (such as biological 
neural networks), some of the neurons may be in quiescent ( "non- firing" ) regime. A more general description of this 
system should incorporate equations for complex amplitudes of oscillators similar to Ref. |12| . 

Our results may have important biological implications. In particular, they suggest that there are fundamental 
reasons for the wide separation of time scales of fast neuronal dynamics ( e.g. membrane potential oscillations) and 
slow synaptic variability. They also present a very simple (possibly, the simplest) model of the synaptic reentry 
reinforcement which is believed to be the foundation of the long-term memory stability |^ . 

We are indebted to M.I.Rabinovich and R. Huerta for fruitful discussions. This work was supported by the U.S. 
Department of Energy grant DE-FG03-96ER14592 and by UC-MEXUS/CONACYT. 
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